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1.  INTRODUCTION 


The  prediction  of  sound  radiation  from  vibrating  bodies  is  a  recurrent 
problem  of  practical  importance  in  engineering  science.  Analytical  solutions 
for  such  problems  are  generally  limited  to  cases  for  which  the  surface  of  the 
object  conforms  to  a  suitable  coordinate  system  such  that  the  wave  equation 
can  be  separated.  For  objects  with  nonstandard  shapes,  an  approach  commonly 
used  is  to  reformulate  the  problem  as  a  boundary  integral  relation,  known  as 
the  Ki rchhof f-Hel mhol tz  integral  theorem  [Pierce,  1981],  which  enables  one  to 
express  the  acoustic  field  at  any  point  exterior  to  the  object  as  a  definite 
integral  over  the  surface  of  the  object.  Terms  in  the  integrands  involve  both 
the  pressure  and  the  normal  acceleration  of  the  surface.  These  two  surface 
quantities,  however,  are  not  independent  and  cannot  be  prescribed  simul¬ 
taneously.  In  most  acoustic  problems,  one  of  the  two  is  specified  and  it  is 
necessary  to  solve  the  integral  equation  for  the  other  unknown  surface  quan¬ 
tity. 

Many  numerical  methods  have  been. developed  in  that  regard,  for  example,  by 
Chen  and  Schweikert  [1963],  Banaugh  and  Goldsmith  [1963],  Chertock  [1964], 
Copley  [1967,  1968],  Schenck  [1968],  Burton  and  Miller  [1971],  Bell  et  al 
[1978,  1979]  for  acoustic  radiation  and  scattering  problems.  The  present 
report  describes  a  variational  formulation  of  the  fluid-structure  interaction, 
which  is  derived  from  the  Ki rchhof f-Hel mhol  tz  integral  theorem.  The  varia¬ 
tional  principle  is  attractive  in  that  it  always  selects,  from  among  various 
adjustable  parameters  of  a  chosen  set  of  approximating  functions,  the  ones 
that  are  optimal.  It  retains  the  accuracy  of  the  boundary  integral 
approaches,  while  offering  the  computational  efficiency  of  approximate 
decoupling  techniques,  such  as  doubly  asymptotic  approximation  developed  by 
Geers  [1971]. 

Although  the  formulation  of  an  integral  equation  as  a  variational  princi¬ 
ple  dated  back  as  early  as  in  1884  [Volterra],  it  was  not  used  for  practical 
calculations  in  wave  diffraction  or  scattering  until  the  late  1940's.  Levine 
and  Schwinger  [1948,  1949]  employed  a  variational  principle  to  study  the  dif¬ 
fraction  of  a  scalar  plane  wave  by  an  aperture  in  an  infinite  plane  screen. 
Levine  [1950],  Bouwkamp  [1954],  and  Sleator  [1960]  further  discussed  varia¬ 
tional  principles  for  acoustic  diffraction  and  scattering  problems. 
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In  general,  derivation  of  a  useful  variational  principle  from  an  integral 
equation  requires  that  the  integration  kernel  be  symmetric,  and  therefore 
self-adjoint.  The  variational  principle  derived  by  Morse  and  Feshbach  [1953] 
was  obtained  by  simply  taking  a  normal  derivative  of  the  Kirchhoff-Helmholtz 
integral  relation  for  the  pressure  at  an  exterior  point,  and  then  letting  the 
exterior  point  approach  the  surface.  The  resulting  principle  features  an 
integration  kernel  that  is  self-adjoint.  However,  one  of  the  resulting 
integrands  becomes  highly  singular,  and  the  integral  is  ambiguous  without 
careful  definition  of  how  the  singular  integral  is  to  be  evaluated.  It  is  not 
apparent  how  their  principle  could  be  numerically  implemented  for  cases  of 
interest. 

In  order  to  circumvent  this  singularity  problem,  the  present  report  util¬ 
izes  a  method  previously  employed  by  Maue  [1949]  and  Stallybrass  [1967].  By 
using  some  mathematical  identities  associated  with  the  free-space  Green's 
functions,  the  integrand  with  non-integrable  singularity  is  recast  into  a  form 
involving  the  tangential  derivative  of  the  surface  pressure.  The  singulari¬ 
ties  contained  in  the  resultant  integrands  are  at  most  of  Cauchy  type 
[Churchill,  1960],  and  therefore  well  behaved  in  the  limit  as  the  external 
point  approaches  the  surface. 

After  we  derive  the  present  variational  principle  for  an  arbitrary  body, 
we  shall  specialize  it  to  situations  where  the  body  is  slender,  such  that 
wetted  surfaces  are  proximate.  The  application  of  the  present  variational 
principle  to  a  circular  disk  with  infinitesimal  thickness  in  rigid  transverse 
vibration  is  demonstrated.  The  numerical  results  agree  remarkably  well  with 
those  obtained  by  Leitner  [1949]. 

2.  KIRCHHOFF-HELMHOLTZ  INTEGRAL  THEOREM 

The  variational  principle  discussed  in  the  present  report  is  based  on  the 
standard  Kirchhoff-Helmholtz  integral  theorem  [Pierce,  1981]  for  the  complex 
pressure  amplitude  p  at  an  external  field  point  jfc  due  to  monochromatic 
excitation  on  a  closed  surface  S.  It  is  an  exact  corollary  of  the  Helmholtz 
wave  equation  and  the  Sommerfeld  radiation  condition  that,  for  *  external  to 
S,  one  has 
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p(x)  -  - 


n(|)*V^p(|)  G(xlf)  dA^ 


4* 


P(?)  n(?)*V^G(xl|)  dA^ 


(1) 


where  £  is  a  source  point  on  S;  the  operator  denotes  the  gradient  at  that 
location,  and  dA^  is  an  element  of  area  surrounding  the  source  point.  The 
quantity  G(£l£)  is  the  free  space  Green's  function  for  the  pressure  at  the 
field  point  t  associated  with  a  point  source  at  £  oscillating  as  e"1ut;  that 
is, 


G(xl|)  =  i  eikR  (2) 

where  k  =  u/c  (with  c  denoting  the  ambient  speed  of  sound),  and  R  is  the 
distance  between  the  source  and  field  points, 


R  ■  lx  -  ?l  (3) 

Our  interest  here  shall  be  focused  on  radiation  problems,  such  that  the 
surface  velocity  is  presumed  a  given  quantity,  while  the  acoustic  pressure  is 
to  be  calculated.  Euler's  equation  of  motion  for  a  fluid  (the  fluid  dynamic 
counterpart  of  Newton's  second  law),  combined  with  continuity  of  the  normal 
component  of  particle  velocity  on  the  surface  leads  to 


n(J)*V^p(|)  =  i  w/n/n(?)  (4) 

where  vn(£)  is  the  complex  amplitude  of  the  (outward)  normal  component  of  the 
surface  velocity. 

Equations  (1)  and  (4)  show  that,  if  the  pressure  and  normal  velocity  on 
the  surface  S  are  known,  then  the  pressure  at  any  external  point  not  on  S  may 
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be  obtained  from  a  definite  integral.  Our  concern  in  this  study  is  the 
evaluation  of  the  pressure  on  S  corresponding  to  the  velocity  distribution 
vn (|) •  Although  we  shall  derive  a  different  interrelationship  between  these 
two  physical  quantities  on  the  surface,  it  is  instructive  to  first  review  the 
derivation  of  a  better  known  relationship  [Schenck,  1968]  that  has  often  been 
used  in  earlier  numerical  studies.  This  relation  is  formally  obtained  by 
letting  t  approach  the  surface  S. 

3.  FIRST  INTEGRAL  EQUATION  FOR  SURFACE  PRESSURE 

In  the  situation  where  t  lies  on  S,  such  that  t  =  £,  both  integrals  in 
Eq.  (1)  have  a  singularity  at  £  =  £  because  R  =  0.  The  singularity  in  G(*l£) 
is  integrable,  but  evaluation  of  the  integral  containing  V£G(*I£)  requires 
careful  consideration.  In  particular,  one  finds  that  the  value  of  the 
integral  may  have  different  values  depending  on  whether  one  regards  t  as 
having  approached  the  surface  from  the  exterior  or  from  the  interior. 

The  analysis  here  of  the  effect  of  the  singularity  is  similar  to  that  of 
Kellogg  [1953],  in  which  the  external  point  *  is  brought  to  a  location  £  on 
the  surface  in  a  limiting  operation.  As  shown  in  Fig.  1,  point  £  is  defined 
to  be  concurrent  with  the  normal  that  intersects  t,  so 


x  =  J  +  e  n(f)  (5) 

where  e  is  the  small  perpendicular  distance  of  t  from  the  surface.  The  region 
S'  is  a  circular  segment  of  S,  centered  at  £,  with  small  radius  a.  The 
remainder  of  the  surface  is  denoted  as  S".  The  required  integral  is  evaluated 
by  taking  the  limit  as  e  0  with  a  fixed,  and  then  taking  the  limit  as  a  -*■  0. 
The  order  in  which  the  two  limits  are  taken  is  important. 

Given  that  e  and  a  are  both  small  compared  with  any  characteristic  dimen¬ 
sions  of  the  surface  S  and  given  that  S  is  smooth  near  the  point  £,  it  is 
appropriate  to  regard  S'  as  having  the  shape  of  an  elliptical  bowl.  The 
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Figure  1.  Integration  when  the  field  point  approaches  the  surface. 
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principal  curvature  radii,  Rj  and  Rn,  of  this  bowl  are  those  of  the  surface 
at  location  £.  In  a  local  coordinate  system  with  origin  at  £,  the  surface  S' 
is  locally  described  by 


z  2  -  x2/(2Rj)  -  y2/(2Rn)  (6) 

(A  convex  surface  would  correspond  to  negative  radii  of  curvature.)  The  unit 
outward  normal  vector  fl(f)  at  a  generic  point  £  on  S'  is  given  approximately 
by 


n(|)  =  ez  +  (x/Rj)ex  +  (y/Rn)ey  (7) 

In  this  local  coordinate  system,  the  external  point  t  is  at  etz,  while  the 
z-coordinate  of  £  is  as  described  by  Eq.  (6)  above.  Consequently,  the  vector 
is  given  by 


R  -  x-|  2  [c  +  (x2/2Ri)  +  (y2/2Rn)]ez  -  xex  -  yey 
and  one  derives 


(8) 


R  n(|)*V^R  s  -€  +  (x2/2Rj)  +  (y2/2Rn)  (g) 

n(?)*V^6(|lx)  =  (1/R'3  +  k2/2R)  [e  -  (x2/2Rj)  -  y2/2Rn]  (io) 

where,  in  the  latter  expression,  the  neglected  terms  are  of  the  order  of  unity 
or  smaller  when  e  and  a  are  both  much  smaller  than  Rj  and  Rjj.  Moreover,  on 

the  right  side  of  Eq.  (10),  it  is  consistent  to  replace  the  distance  R,  wher¬ 
ever  it  appears,  by 


R  -  (e  +  r  ) 


.2a/2. 


2  2  2 
r  =  x  +  y 


(11) 
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and,  for  the  differential  of  area,  to  set 


dA^  =  rdr  d^;  x  =  r  cos^;  y  =  r  sin^  (12) 

where,  for  the  surface  S',  the  r-integration  extends  from  0  to  a,  and  the 
^-integration  extends  from  0  to  2t. 

Since  all  points  on  S"  are  at  a  finite  distance  from  t,  the  factor 
fl(£)*V£G(*l|)  has  no  singularities  on  S"  when  e  +  0.  Moreover,  if  one  lets 
e  0  first  and  considers  r  to  be  small  but  nonzero,  then 


n(?)*V^G(xl|)dA^  -*•  -  (1/2)  (Rj*cos2^  +  R'Jsin2^) (1/r)  rdr  d<j>  (13) 

The  1/r  singularity  here  is  cancelled  by  the  factor  in  the  area  differential 
factor  rdr.  The  integral  in  this  limit  is  exactly  the  same  as  if  one  set  t  to 
£  at  the  outset,  and  then  did  the  integral  over  the  surface  S.  Hence,  it  is 
mathematically  meaningful  to  integrate  over  the  surface  S"  by  taking  the 
double  limit  such  that  first  e+0,  then  0.  We  therefore  turn  our  attention 
to  the  contribution  of  S'. 

Let  F(£)  be  any  continuous  scalar  function  on  S.  Then  Eqs.  (10)  and  (12) 
lead  to 


F(|)  n(f)*V£G(xl?)  dA 

J  J  *  * 

S' 

-  F(J)  (e/R3)  rdr  d ^ 

J  d 

S' 


F(J)  [A(e,r,^)  -  (r/R)2B(^)] R_1rdr  d0 


(14) 
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where  in  the  latter  term  we  abbreviate 


(15a) 

(15b) 


A(e,rf0)  =  (l/2)k2e  -  (l/4)k2r2(Rj1cos2^  +  Rjjsin2^) 

B(^)  2  (l/2)(Rj1cos2^  +  Rjjsin2^) 

Functions  A  and  B  are  both  finite;  A  actually  vanishes  in  the  limit  as  both  e 

and  a  go  to  zero.  Also,  we  note  that  the  quantity  r/R  is  always  less  than  1, 

regardless  of  the  value  of  e.  Consequently,  the  second  term  in  Eq.(14)  is  of 
the  order  of  magnitude  of  a\  it  therefore  vanishes  in  the  limit  a  -*•  0. 

We  now  proceed  to  evaluate  the  first  term  in  Eq.  (14).  With  e  small, 
but  positive,  an  appropriate  change  of  integration  variable  is  to  the  polar 
angle  9,  defined  such  that 

r  =  e  tan0;  R  =  e  sec0;  dr  =  e  sec 20  d0;  (e/R3)  r  dr  =  sin0  d 9  (16) 

The  integration  limits  on  9  are  0  and  tan_l(a/e) .  In  the  limit  as  e+0  with  a 

fixed,  this  upper  limit  approaches  t/2.  The  substitution  of  Eq.  (16)  into 
Eq.  (14)  makes  the  integration  trivial,  the  overall  result  being  simply 
2*‘F(^‘)«  Thus,  we  find  that,  as  regards  the  integral  over  the  entire  surface, 
when  e  goes  from  a  finite  positive  value  to  0, 

{  * 

l™  o  F(?)  n(?)*V,G(xl|)  dA, 

J  J  *  * 

S 

=  2t  F($-)  +  F(|)  n(|)*V.G(J| J)  dA*  (17) 

J  J  *  > 

S 

The  result  of  applying  Eq.  (17)  to  the  Kirchhoff-Helmholtz  integral  equa¬ 
tion  (1),  with  the  pressure  boundary  condition  given  by  Eq.  (4),  is 
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p  (?) =  - 


r  f 


S 


iu/>vn(|)  G(Jl|)  dA^ 


S 


P(f)  n(?)*V$G(?l|)  dA^ 


(18) 


The  above  integral  equation  has  been  used  often  for  numerical  analysis  of 
surface  pressure,  but  the  numerical  implementations  have  not  led  to  especially 
accurate  results  for  many  cases  of  practical  interest.  One  of  the  intrinsic 
sources  of  difficulty  [Rogers,  1973]  is  that,  for  certain  discrete  character¬ 
istic  frequencies,  the  solution  is  not  unique,  because  the  homogeneous 
integral  equation  has  eigensolutions  (the  multiplicative  constant  being 
arbitrary)  at  these  frequencies.  Another  complication  is  that  the  integration 
kernel,  li(^)*V^G(^l £) ,  although  integrable,  is  singular,  so  one  is  confronted 
with  a  singular  integral  equation.  Also,  the  kernel  is  not  symmetric  in  the 
interchange  of  £  and  £,  so  the  linear  operator  associated  with  this  integral 
equation  is  not  self-adjoint. 


4.  SECOND  INTEGRAL  EQUATION  FOR  SURFACE  PRESSURE 


A  strong  case  can  be  made  that  a  variational  principle  is  often  a  useful 
intermediate  step  for  developing  robust  numerical  methods  for  solving  an 
integral  equation.  However,  the  integral  equation  (18)  does  not  lead  in  a 
natural  manner  to  a  variational  principle,  because  the  integration  kernel  is 
not  self-adjoint.  Consequently,  we  have  made  use  of  a  second  integral  equa¬ 
tion  (strictly  speaking,  a  differential-integral  equation)  that  applies  to  the 
same  problem  [Meyer  et.  al.,  1978].  Although  this  equation  is  less  well 
known,  it  nevertheless  has  the  standard  Kirchhoff-Helmholtz  integral  theorem 
as  its  foundation.  We  begin  by  using  Eq.(l)  to  form  the  normal  derivative  of 
the  pressure  on  the  surface.  Because  of  the  singularity  of  the  integral  terms 
for  a  field  point  at  the  surface,  we  shall  temporarily  keep  t  off  S  by  making 
use  of  Eq.  (5).  Thus,  we  form 
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'  fe  n(J)*VxG(xl|)  dA. 

J  J  ^ 


S 


fell  P(?)  Cn(J)*Vx]  [n(|).?e]  G(Slf)  dAe 


(19) 


S 


where  Vx  denotes  the  gradient  at 

Note  that  the  integrand  in  the  second  term  of  Eq.  (19),  in  the  limit  as 
1S  highly  singular,  so  the  appropriate  expression  for  numerical  eval¬ 
uation  of  this  term  in  the  limit  of  e  +0  has  to  be  developed  with  some  care. 
To  this  purpose,  we  follow  an  approach  similar  to  that  previously  developed  by 
Maue  [1949]  and  Stallybrass  [1967],  We  first  introduce  the  vector  identity 


(a  x  8)  •  (c  x  3)  «  (a  •  c)  (S  •  3)  -  (a  •  3)  (S  •  c) 


(20) 


so  that 


(21) 


where  Sn-  and  are  unit  vectors  (possibly  the  same)  in  any  cartesian 
coordinate  system  and  G  is  the  free-space  Green's  function  (2).  Since  Gtxlf) 
is  a  function  only  of  R  =  V£G(*l£)=-VxG(*l£) .  Thi  s  allows  us  to  replace 
the  third  term  on  the  right  side  of  the  above  equation  by  an  analogous 
expression,  in  which  the  subscripts  i  and  j  are  interchanged.  After  making 
such  a  change,  we  multiply  each  term  by  nj(£)nj(£),  and  then  sum  over  i  and  j 
from  1  to  3.  A  rearrangement  of  terms  leads  to 
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(22) 


[n(?)-Vx]  [n(?)*Ve]G  =  [n(J)*n(?)]  V^G 


-  [n(J)  x  Vx]*[n(|)"x  V^]G 

Also,  because  G(*l|)  has  the  property  that  V^G(jtl^)  is  the  negative  of 
^xG(*l t) i  an^  because  it  satisfies  the  scalar  Helmholtz  equation,  we  can  make 
the  substitution 


Vx*7eG(xl|)  =  -  V2G(xl|)  =  k2G(xl|) 

These  two  latter  relations  allow  Eq.  (19)  to  be  recast  in  the  form 


(23) 


n(J)*VYp(x)  =  -  ^ 


iw/w,/!)  n(J)*VxG(xl|)  dA 


S  . 


P(|)G(xl|)n(|)  dA 


-  j;  [S(f)  x  7  ]• 


P(|)[n(|)  x  V£]  G(xlf)  dA 


(24) 


The  second  derivative  of  the  Green's  function  will  present  a  problem  when  e+0, 
so  we  integrate  the  last  term  by  parts,  such  that 
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P(£) [n(|)  x  Ve]  G(xlJ)  dA^ 


S 

=  [n(?)  x  VJ  {ptctG(xt$T}  dA* 

J  J  *  * 

S 

^  * 

[n(|)  x  V>p(|)]  G(xl|)  dA*  (25) 

J  J  * 

S 

The  first  integral  on  the  right  side  vanishes.  In  order  to  demonstrate  this, 
we  recall  a  version  of  Stokes'  theorem  [Less,  1950].  Let  F  be  any  scalar 
function  that  is  continuous  and  that  is  differentiable  with  respect  to  tangen¬ 
tial  coordinates  in  S",  and  let  C  denote  the  curve  bounding  S".  Then 


n(|)  x  V/(|)  dA, 

J  J  *  * 

S" 


(26) 


where  the  line  element  dt  progresses  along  C  in  the  counter-clockwise  sense 
when  the  unit  normal  fl  points  toward  the  observer.  In  the  present  situation, 
F  represents  the  factor  contained  within  braces  in  the  first  integral  in 
Eq.  (25).  Since  e*0,  this  factor  F  is  finite.  Also,  since  S  is  a  closed 
surface,  one  can  regard  the  open  surface  S"  as  being  all  of  S  except  for  a 
small  patch  surrounding  any  arbitrarily  chosen  point  on  S;  in  the  limit  as  the 
dimensions  of  this  patch  shrink  to  zero,  the  curve  C  shrinks  to  zero  length. 
It  follows  that  the  right  side  of  Eq.  (26)  vanishes  when  S"  becomes  S.  As  a 
result,  Eq.  (24)  becomes 


-  15  - 


n(J)-?j*(x)  = 


h 


vn(J)  n(J)«?xG(xl  \)  dA 


Tr  «©• 


n(?)  G(xl|)  p(|)  dA 


S 


e 


+  b  £*©  *  V 


G(xl|)  n(|)  x  V$p(|)  dA^ 


(27) 


The  desired  differential -integral  equation  emerges  after  one  takes  the 
limit  of  Eq.  (27)  as  e  +  0;  doing  so  yields 


1w/>un(^  =  h 


n(|)  G(J||)  p(|)  dA 


e 


+  G®  x  V^].jjG(Jl|)  n(^)  x  Vep(?)  dAe 


(28) 


where  we  abbreviate 


u„6>  ■  vn(J>  +  e'™  0  4?  v„(?)  dA 


K®  +  fcPR  *„(?>  S6>*¥;Jr 


ikR 


dA, 


(29a) 


(29b) 


In  the  latter  version,  the  symbol  PR  implies  principal  value  and  R  is  under¬ 
stood  to  be  l£-fl.  The  equivalence  of  these  two  versions  of  Eq.  (29)  can  be 
demonstrated  in  a  manner  similar  to  that  used  in  deriving  Eq.  (17).  The 
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second  version,  however,  has  limited  use;  it  should  not  be  used,  for  example, 
if  the  vibrating  body  is  an  unbaffled  plate  (curved  or  not  curved)  with 
thickness  idealized  as  being  infinitesimal. 

5.  DERIVATION  OF  THE  VARIATIONAL  PRINCIPLE 

We  shall  derive  the  variational  principle  from  Eq.  (28).  The  first  step 
is  to  multiply  each  term  of  that  equation  by  a  virtual  increment  <5p(£)  in  the 
actual  pressure  amplitude  at  the  surface  point  £.  Each  term  is  then 
integrated  over  the  surface  S,  with  the  parameter  £  distinguishing  the  points 
on  the  surface  during  this  integration.  The  result  is  then  that 


i  up 


*P(J)  Un(J)  dA< 


j c 

4  IT 


5p(J)  [n(J)*n(?)]  p(|)  G(J||)  dA^dA^. 


S  S 


4t 


<?P(£)[n(£)  x  JjG(Jl|)[n(?)  x  V^]p(|)  dAj  dA^.  (30) 


An  analysis  (involving  Stokes'  theorem)  similar  to  that  described  in  Eqs.  (24) 
and  (25)  applies  to  the  last  term  in  Eq.  (30)  with  the  inner  integral  (en¬ 
closed  in  large  braces)  regarded  as  a  function  F(£).  Such  a  procedure  enables 
one,  in  effect,  to  integrate  by  parts,  transferring  (with  a  change  in  sign) 
the  operator  rt(£)  x  to  the  factor  <5p(£).  The  result  is  that 
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1  up 


tfp(J)  un($)  dA< 


4t 


*P(J)  [n(J)*n(?)]  p(|)  G(Jlf)  dA^ 


S  S 


1_ 

4t 


[n(J)  x  V^5p(J)]*[n(f)  x  V^p(|)]G(Jl|)  dA^dA^. 


(31) 


S  S 


We  are  seeking  the  pressure  amplitude  p  at  points  on  the  surface  S,  given 
the  normal  velocity  amplitude  distribution  vn(f).  This  means  that  vn(|) 
should  be  held  constant  when  p  is  varied.  Symmetry  of  the  Green's  function 
and  the  usual  rules  of  the  calculus  of  variables  makes  it  possible  to  replace, 
within  the  integrands,  variational  factors  of  the  generic  form  F(|)<5F(£)  by 
alternative  variational  factors  (l/2)$[F(t)F(£)] .  Also,  because  the 
integration  limits  do  not  change  during  the  variation,  the  integral  of  the 
variation  is  the  variation  of  the  integral.  Similarly,  a  sum  of  variations 
can  be  replaced  by  the  variation  of  the  sum.  All  this  allows  us  to  recast 
Eq.  (31)  as  the  variational  principle 


<5J  [p]  =  0 

where  the  functional  J[p]  can  be  identified  as 


(32) 
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J[p]  =  -  iw/7  p(J)  4TUn(J)  dA 

JJ  ** 

s 

2  f  f  f  f 

+  Y  [n (J) •"(?)]  P(J)  P(?)  G(Ji?)  dA.dA 

J  J  J  J  V  j 

s  s 

(  (  C  f 

-  \  [n(J)  x  Vj)(J)]*[n(|)  x  V.p(|)]G(J||)  dA.dA  (33) 

J  J  J  J  J  *  »  J 

S  S 

The  implication  of  Eq.  (32)  is  that  the  functional  J  is  stationary  to 
small  changes  <5p(£)  in  the  pressure  distribution  p(£).  One  cannot,  however, 
in  general  state  that  J  has  an  extreme  value  (maximum  or  minimum)  when  the 
trial  function  p(£)  is  equal  to  the  true  complex  pressure  distribution.  If 
the  trial  function  is  taken  to  be  the  actual  function  plus  ef(£),  where  e  is 
small  and  f(£)  is  a  fixed  admissible  function,  then  J [p]  can  be  regarded  as  a 
function  of  e,  and  dJ/de  must  be  zero  at  e=0.  The  sign  of  either  the  real  or 
imaginary  part  of  the  second  derivative  d2J/de2,  however,  depends  on  the 
choice  of  this  fixed  admissible  funtion  f(£).  In  general,  the  class  of 
admissible  trial  functions  is  restricted  to  functions  that  are  continuous  over 
the  surface  and  piece-wise  differentiable  with  respect  to  displacements  on  the 
surface,  such  that  fi(£)  x  V^p(£)  exists  almost  everywhere  and  is  square 
integrable. 

One  possible  trial  function  is  simply  a  position-independent  factor  times 
the  actual  function.  The  stationarity  condition  is  satisfied  when  this  factor 
is  unity,  but  differentiation  with  respect  to  this  factor  introduces  addi¬ 
tional  factors  of  1  and  2,  respectively,  for  those  terms  in  Eq.  (33)  which  are 
linear  and  bilinear  in  the  surface  pressure.  Consequently,  after  differen¬ 
tiating  the  functional  with  respect  to  the  factor  and  then  setting  the  factor 
equal  to  1,  one  derives 
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107  p 


+  k* 


P(J)  4jrUn(J)  dA( 


[n(J)*n(|)]  p(J)  p(|)  G(J|J)  dA^dA^. 


S  S 


-  J  Cn(J)  x  V^p(J)] •[«(!)  x  7^p(|)]G(Jl|)  dA^dA  (34) 
s  s 

where  the  inserted  function  p(£)  is  the  actual  complex  pressure  on  the  sur¬ 
face.  (Note  that  the  foregoing  is  merely  the  integral  over  dA^  of  Eq.  (28).) 
Substitution  of  Eq.  (34)  into  Eq.  (33)  reveals  that  the  stationary  value  of 
the  functional  J  is  numerically  equal  to 


J[ptrue]  =  ‘  2  ptrue&  4xUn®  dA< 

J  J  * 

S 


(35) 


Consequently,  if  a  trial  function  for  p(£),  that  is  correct  to  first  order,  is 
inserted  into  the  functional  J[p]  stated  in  Eq.  (33),  one  obtains  an  estimate, 
accurate  to  second  order,  of  the  surface  integral  of  Ptrue(t)Mt)  •  In  some 
cases,  the  latter  may  be  of  principal  interest  and  may  have  an  important 
physical  identification;  an  accurate  method  for  estimating  its  value  would 
then  be  of  great  use. 


6.  INTEGRALS  INVOLVING  SURFACE  VELOCITY 

A  few  guidelines  may  be  stated  at  the  outset  concerning  the  improper  inte¬ 
gral  that  causes  the  distinction  between  Un  and  vn  appearing  in  Eqs.  (29).  We 
rewrite  those  relations  as 
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(36a) 


un©-*n© 


1  im 
e  0 


1_ 

4r 


vn(?)  n(J)«VxG(xl|)  dA£ 
S 


-  5  »n®  +  fc  *fk<t> 


Ht 

R  dR 


ikR 


dA, 


(36b) 


For  thin  platelike  structures  that  are  unbaffled,  the  surface  S  consists 
of  two  contiguous  sheets  with  infinitesimal  separation  between  them.  If  these 
sheets  are  labelled  by  the  subscripts  I  and  II,  then  adjacent  points  on  the 
front  and  back  of  the  surface  can  be  labelled  and  <£jj.  The  geometry 
requires 


n(|j)  =  -n(?n) 

so,  if  the  surfaces  are  vibrating  together  as  a  unit,  one  must  have 


(37) 


vn^i)  "  'vn^II^  (38) 

Insofar  as  the  integral  on  the  right  side  of  Eq.  (36a)  is  concerned,  the 
separation  between  and  is  much  less  than  e  while  the  limit  is  being 
taken.  Consequently  Rj  and  Rji  are  the  same  for  all  finite  e;  the  opposite 
signs  in  Eq.(38)  therefore  require  the  integrations  over  the  two  sides  of  S  to 
be  equal  but  opposite  for  all  finite  e;  these  integrals  cancel  and  one  is  left 
with  the  simple  result 


n  =  vn  (39) 

for  thin  unbaffled  plate-like  bodies  (alternately  referred  to  as  laminas). 
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For  less  specialized  circumstances,  a  fail-safe  method  for  evaluating  the 
integral  in  Eq.  (36a),  without  the  complications  of  taking  a  limit  or  of 
separately  distinguishing  circumstances  when  portions  of  the  body  are  laminas, 
can  be  developed  from  the  identity 


e  ?  0  "(tJ’VR'1)  dA  =  0  (40) 

J  J  * 

s 

which  follows  from  (i)  the  fact  that  Vx  applied  to  1/R  is  equivalent  to  -7^ 
applied  to  1/R,  (ii)  Gauss's  theorem,  and  (ill)  the  fact  that  1/R  satisfies 
Laplace's  equation.  Consequently,  one  can  write 


4tK  (?>-%©] 


1  im 
e+0 


Cvn(t)n(J) 

-J  J 
S 


v„(f)  n©*VxG(xl|)  dA^ 


S 


©n(?)3*  ^ 


d_ 

dR 


_  i  kR 
e _ 

R 


JJ 

s 


ikR 


1 

R 


dA, 


(41) 


Here  the  singularities  in  the  integrands  of  the  integrals  in  the  second 
version  are  at  most  of  order  1/R  (which  is  integrate) ,  so  one  need  not  be 
concerned  with  explicitly  taking  the  limit  as  e+0. 


An  example  for  which  the  above  representation  might  be  computationally 
useful  is  that  of  a  finite  length  circular  cylinder  vibrating  as  a  rigid  body 
parallel  to  its  axis.  If  the  length  L  of  the  cylinder  goes  to  zero,  then  one 
has  a  circular  disk  (discussed  in  greater  depth  in  the  next  section)  in 
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transverse  vibration,  and  Eq.  (39)  should  apply,  Un  =  vn.  This  limit,  how¬ 
ever,  does  not  emerge  very  easily  from  Eq.  (36b),  but  it  does  pop  out  of 
Eq.  (41).  In  the  limit  of  L+0,  the  area  integrals  in  (41)  are  over  only  the 
top  (z  =  L/2)  and  bottom  surfaces  (z  =  -  L/2)  of  the  cylinder.  If  ^  is  a 

point  on  the  top  surface,  then  fl(£)  is  in  the  direction  of  the  unit  vector  tz, 
while  rt(£)  may  be  in  either  the  $z  or  -tz  directions.  With  L-K),  one  has  in 
any  case  that  is  perpendicular  to  both  fi(£)  and  rt(£),  so  both  the  first 
and  second  terms  in  Eq.  (41)  vanish  trivially  in  this  limit.  One  cannot  draw 
such  a  conclusion  for  the  integral  in  Eq.  (36b),  because  the  principal  value 
applies  only  for  the  integral  over  the  top  surface  (given  that  £  is  on  the  top 
surface).  For  the  integral  over  the  bottom  end,  one  must  in  general  evaluate 
that  integral  for  finite  L,  then  take  the  limit  as  L+0.  One  would  get  a 
different  (and  incorrect)  answer  if  one  jumped  inside  the  integral  and  took 

the  limit  as  L+O  before  evaluating  the  limit.  Such  problems  do  not  arise, 

however,  for  the  integrals  in  the  latter  version  of  Eq.  (41) ,  because  the 
integrands  are  sufficiently  well-behaved  and  non-singular  that  the  order  of 
taking  the  limit  and  carrying  out  the  integrations  can  be  freely  interchanged. 

7.  SPECIAL  CASE  OF  A  THIN  DISK 

Because  there  are  no  gradients  of  the  Green's  function  within  the  inte¬ 
grands  of  the  second  and  third  terms  of  the  functional  J[p]  in  Eq.  (33),  one 
need  not  characterize  those  integrals  as  being  Cauchy  principal  values.  Con¬ 
sequently,  the  formalism  is  easily  adapted  to  a  slender  body,  for  which  one 
region  of  the  surface  S  is  infinitesimally  close  to  a  different  region  of  S. 

For  example,  consider  the  disk  in  Fig.  2,  which  is  shown  on  edge.  Let  us 
denote  variables  on  the  upper  and  lower  surfaces  by  a  subscript  plus  or  minus 
respectively.  The  disk  is  assumed  to  be  vibrating  as  a  rigid  body  normal  to 
its  faces,  so  the  analysis  given  in  the  preceding  section  applies  and  one  has 

un  =  vn*  Also,  as  stated  in  Eq.  (38),  the  surface  velocities  vn(|+)  and 

vn(t~)  are  180*  out-of-phase,  that  is,  when  the  upper  surface  moves  inward, 
the  lower  surface  moves  outward,  and  vice  versa.  It  follows  that  the  actual 
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field 
point  x 


point  r+ 


Figure  2.  Field  and  surface  points  for  a  thin  disk. 
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acoustic  pressure  must  have  comparable  behavior,  and  must  be  antisymmetric 
with  respect  to  reflection  through  the  plane  on  which  the  disk  nominally  lies. 
We  choose  to  restrict  our  class  of  admissible  trial  functions  for  p  to  include 
only  functions  that  have  this  property.  Thus  we  have 


-  -vn(?+).  P(?J  =  -P(?+)  (42) 

%(?>„(?_)  ■  »„(?+)P(?+)  (43) 

Also,  since  the  unit  normal  vectors  on  opposite  sides  of  the  disk  are 
oppositely  directed,  we  have 

"(£>(?->  ■  »(tjp(£js  "(?.)  X  TjP (tj  =  n(J+)  X  V?p(?+)  (44) 

These  relations,  of  course,  also  hold  if  the  dummy  integration  variable  £  is 
replaced  by  the  dummy  integration  variable  £. 


Equations  (42-44)  substantially  simplify  the  functional  J[p].  When  each 
of  the  integrals  in  Eq.  (33)  is  decomposed  into  the  contributions  of  the  upper 
and  lower  surfaces,  we  find  that  J[p]  may  be  expressed  in  terms  of  integrals 
extending  over  the  upper  surface  only;  specifically, 


J[p] 


Sriup 


(?)p(?)  dA^  +  2k2 

J  4 
+ 


(  f 


J 


P(J)P(?)G(Jl|)  dA^dA^. 


&zx  (J) ] •  Cezx  V^p(|)]G(Jl|)  dA^dA 


S+  S+ 


where  £  and  £  are  dimensionless  variables  and  can 


(45) 
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8.  ANALYSIS  OF  A  RIGID  DISK  WITH  DISTRIBUTED  BASIS  FUNCTIONS 

In  this  section  we  shall  represent  the  pressure  distribution  on  the 
surface  of  an  unbaffled  rigid  disk  in  a  series  of  assumed  basis  functions, 
which  are  each  in  general  nonzero  at  every  point  on  the  disk.  It  is  conven¬ 
ient  to  begin  by  nondimensionalizing  the  functional  J  in  Eq.  (45).  For  this 
purpose  we  scale  distance  by  the  radius  a  of  the  disk,  and  scale  pressure  by 
/jcv0,  where  p  and  c  are  the  ambient  density  and  sound  speed,  and  v0  is  the 
amplitude  of  the  surface  velocity,  such  that  on  the  plus  side  of  the  disk,  the 
normal  velocity's  complex  amplitude  is  given  by 


vn(^  =  vo  (46) 

Let  an  overcarat  denote  a  nondimensional  quantity.  We  then  find  that  Eq.  (45) 
becomes 


J[p] 


4/c2v^a 


where 


-  2xika 


P(£)  dA.  + 


A  A  A  A  A  A  A  A  a 

P(£)P({)G(£|£)  dA^dA^. 


S+  S+ 


1 

2 


A  A 
S+  S+ 


[ezx  V^p(J)].[e2x  V^p(e)]G(Jie)  dA^ 


(46) 


P  = 


/>cv0* 


V  =  a  V;  G  =  a  G  = 


;ikaR 

A  f 

R 


dA  =  -j  dA 
a^ 


(47) 
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Because  of  the  axial  symmetry  of  the  surface  velocity,  the  pressure 
amplitude  on  the  disk  must  only  be  a  function  of  the  radial  distance  from  the 
center  of  the  disk.,  Polar  coordinates  for  the  disk  are  depicted  in  Fig.  3, 
where  ra  and  sa  are  defined  as  the  radial  distances  to  the  points  £  and  £, 
respectively.  Thus  r  and  s  are  dimensionless  variables  that  may  be  regarded 
as  the  possible  arguments  (which  replace  £  and  £)  of  the  scaled  dimensionless 
pressure.  Correspondingly,  we  find  that 


P(?)  3  P(r);  p(?)  =  p(s); 

P'(s)es; 

(48) 

where  a  prime  indicates  a  derivative  with  respect  to  the  argument,  and  £r  and 
£s  are  the  radial  unit  vectors  shown  in  Fig.  3..  The  nondimensional  distance  R 
is  found  from  the  law  of  cosines  to  depend  on  the  relative  polar  angle  9 
(equal  to  9r-9s) ,  according  to  the  well  known  relation 

R  =  (r2  +  s2  -  2rs  cos 9)^^  (49) 

The  scaled  elements  of  area  in  this  axi symmetric  situation  may  be  taken  as 


=  P'(r)er;  V^p(£) 


AAA 


AAA 

G(fl«)  * 


1  kaR 


dA^.  =  rd0pdr  ,  dA^  =  sd0$ds  (50) 

Since  the  integrand  for  the  double  integrations  over  area  depends  only  on  the 
relative  angle  9,  a  simple  transformation  to  9  as  a  variable  of  integration, 
replacing,  say,  0r,  enables  one  of  the  angular  integrations  to  be  done 
trivially,  yielding  a  factor  of  2t.  This  leads  to  the  following  simple  form 
for  the  scaled  functional  that  appears  in  Eq.  (46) 


-  27  - 


point  £ 


1 


1  1 


A 


-2xi ka 


p(r)r  dr  + 


(ka)2 

2 


0 


p(r)p(s)C^(rl s)  drds 
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2t 

p 

* 
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A 

JkaR 

dd 

R 


2x 


0 


JkaR 
e _ _ 

A 

R 


cos0  dd 


(51) 


(51a) 


(52b) 


Note  that  the  factor  2x  occurs  in  Eq.  (51)  because  the  integrands  in  Eq.  (46) 
are  dependent  on  the  polar  angles  9$  and  9 ^  only  through  their  difference, 
which  is  here  denoted  as  9.  The  differential  factors  Ci(rls)  ds  dr  and 
C2(rls)  ds  dr  can  be  interpreted  as  describing  the  effect  of  an  annulus  of 
radius  s  and  width  ds  on  another  annulus  of  radius  r  and  width  dr.  We  shall 
refer  to  these  functions  Cj  and  C2  as  the  integrated  Green's  functions. 


The  basis  functions  selected  to  represent  the  pressure  distribution  on  the 
surface  of  the  disk  must  satisfy  all  requirements  that  are  imposed  by  the 
manner  of  derivation  of  the  variational  principle  on  admissible  trial  func¬ 
tions.  In  particular,  these  must  be  continuous  and  differentiable,  so  that 
derivatives  such  as  p'(r)  exist.  The  derivation  exploited  the  fact  that  the 
pressure  at  corresponding  locations  on  the  upper  and  lower  surface  differ  only 
in  sign.  However,  because  the  thickness  of  the  disk  is  infinitesimal ,  the 
essure  on  the  two  surfaces  at  the  edge  must  be  equal.  Both  conditions  can 
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(53) 


only  be  satisfied  if  the  pressure  vanishes  at  the  edge,  so 

p(r)  =0  at  r  =  1 

Since  the  admissible  trial  functions  must  be  continuous  over  the  disk's 
surface,  each  of  the  basis  functions  must  satisfy  the  above  condition.  Thus 
the  general  trial  function  formed  as  a  linear  combination  of  N  basis  functions 
can  be  written 


P(r)  «II  Pn*n(r);  *n(l)  «  0  (54) 

Note  that  the  condition  p1  =  0  at  r  =  0,  which  results  from  axi symmetry,  is  a 
natural  boundary  condition  that  will  emerge  from  the  variational  principle  for 
the  exact  solution,  but  it  is  not  a  requirement  that  must  be  imposed  at  the 
outset  on  trial  functions  and  on  variations.  It  is  not  necessary  that  the 
basis  functions  $n(r)  satisfy  this  condition. 

When  we  substitute  Eq.  (54)  and  a  similar  expression  for  p(s)  into 
Eq.  (51),  we  find  that  J  is  a  quadratic  polynomial  in  the  complex  coefficients 
Pn: 
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(55) 
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(56) 


(57) 


The  task  of  evaluating  the  coefficients  ajn  and  bn  shall  be  addressed  in 
the  next  section.  Once  they  are  known,  the  variation  of  J[p]  is  obtained  from 
virtual  increments  in  each  of  the  pressure  coefficients.  Those  increments  are 
arbitrary,  so  <?3  =  0  for  all  admissible  variations  in  the  trial  function  (54) 
requires  that 


gj 

gp-  =  0  for  n  =  1,2, ...,N  (58) 

n 

In  view  of  the  form  of  J  indicated  by  Eq.  (55),  we  conclude  that  the  coef¬ 
ficients  are  governed  by 


[a]{P}  =  {b}  (59) 

where  the  elements  of  the  square  array  [a],  and  of  the  vectors  {b}  and  {P}  are 
the  amn,  bn,  and  Pn,  respectively.  After  the  set  of  linear  equations  repre¬ 
sented  by  Eq.  (59)  has  been  solved  for  the  coefficients  Pn,  it  is  a  simple 
matter  to  recreate  the  spatial  distribution  p(r)  from  Eq.  (54). 
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9.  EVALUATION  OF  THE  INTEGRATED  GREEN'S  FUNCTIONS 


A  singularity  in  each  of  the  integrated  Green's  functions,  given  by 
Eqs.  (52),  arises  when  the  source  radial  coordinate  s  approaches  the  field 
radial  coordinate  r.  In  such  a  circumstance,  the  separation  distance  6,  which 
appears  in  the  denominator  of  the  respective  integrands,  vanishes  for  0=0. 
Our  approach  to  evaluating  the  Cj(rls)  will  separate  out  their  singular  parts 
as  elliptic  integrals. 

The  analysis  and  the  writing  of  the  arguments  of  the  elliptic  integrals 
are  facilitated  if  one  introduces  the  abbreviations 


m  =  — a  =  ka(r+s) 
(r+s)2 

so  that  Eqs.  (52)  become 


(60) 


Cj(rls)  =  j  m(r+s) 


0 


C2(rls)  =  j  m(r+s) 


r 

'  iaD 

-q —  cos0  dd 
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where 


(61a) 


(61b) 


D 
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r+s 


m  cos2(0/2)j 


1/2 


(62) 


The  factor  cos0  in  the  integrand  of  the  expression  (61a)  for  C2(rls)  is 
eliminated  with  the  aid  of  a  trigonometric  identity 
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cosd  =  2  cos2(0/2)  -  1  =  -  -  D2 

mm 

(63) 

As  a  result,  the  integrated  Green's  functions  may  be  obtained 
integrals,  Gj  and  G2,  according  to 

from  two  basic 

C^rls)  =  m(r+s)G1(a,m) 

(64a) 

C2(rls)  =  (r+s)[(2-m)G1(a,m)  -  2G2(<z,m)] 

(64b) 

where 


Gj (a,m) 


G2(a,m) 


x 


0 


(65a) 


(65b) 


Both  of  the  above  expressions  reduce  to  elliptic  integrals  when  a  is  set 
to  zero.  To  highlight  the  resemblance  of  the  functions  Gj(a,m)  to  elliptic 
integrals,  Euler's  identity  and  the  trigonometric  half-angle  formulas  are 
combined.  The  exponential  factor  then  becomes 


e1aD  =  [l  -  2  sin2[f^]]  +  2i  sin  [|^j 
=  1  +  21  sin  exp 

We  now  introduce  a  change  of  variable  for  the 


cos  (r) 

polar  angle  9, 


(66a) 

(66b) 
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(67) 


so  that  the  distance  parameter  D  becomes 


D 


m  si 


1nV) 


1/2 


Substitution  of  Eqs.  (66) -(68)  converts  Eqs.  (65)  to 


(68) 


t/2 


Gj(a,m) 


=  K(m)  +  21 


5  s1n  (r)  ^ppf5)  w 


(69a) 


G2(a,m)  =  E(m)  +  2i 


t/2 

f 

D  sin 


0 


(69b) 


where  K(m)  and  E(m)  are  the  complete  elliptic  integrals  [Milne-Thomson,  1972] 
of  the  first  and  second  kind,  respectively, 


K(m) 


r/2 

d £ _ 

(1-m  sin2^)1//2 

0 


(70a) 


E(m) 


t/2 

(1-m  sin2^)1/,2d^ 

4 

0 


(70b) 


The  singularity  associated  with  the  field  point  approaching  the  source 
point  corresponds  now  to  m  =  1  and  <j>  =  t/2,  in  which  case  D  =  0.  Both 
integrands  in  Eqs.  (69)  have  a  finite  limit  as  D+0,  so  their  evaluation 
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involves  a  straightforward  application  of  conventional  numerical  methods.  We 
use  Simpson's  rule,  based  on  subdividing  the  interval  0<^<x/2  sufficiently 
finely  to  resolve  the  oscillations  in  the  sinusoidal  functions  forming  the 
integrand.  Corresponding  values  for  the  elliptic  integrals  were  obtained  from 
polynomial  approximations  [Milne-Thomson,  1972].  Such  are  conveniently 
expressed  in  terms  of  the  complementary  parameter 


nij  s  1  -  m 


(71) 


Then,  to  an  absolute  error  that  does  not  exceed  2xl0-8,  we  have 


4 


(72a) 


j=0 


4 


(72b) 


j-1 


where  the  numerical  coefficients  are 
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aQ  =  (1/2) ln(16) 

=  1.38629436112 
=  0.09666344259 
a2  =  0.03590092383 
a3  =  0.03742563713 
a4  =  0.01451196212 

cx  =  0.44325141463 
c2  =  0.06260601220 
c3  =  0.04757383546 
c4  =  0.01736506451 


b0  =  0.5 

b1  =  0.12498593597 
b2  =  0.06880248576 
b3  =  0.03328355346 
b4  =  0.00441787012 

d1  =  0.24998368310 
d2  =  0.09200180037 
d3  =  0.04069697526 
d4  =  0.00526449639 


(73) 


The  logarithmic  singularity  in  K(m)  as  m+1  causes  both  of  the  integrated 
Green's  functions  Cj(rls)  to  become  singular  as  s+r.  Although  the  manner  of 
derivation  of  Eq.  (56),  which  defines  the  matrix  coefficient  ajn,  indicates 
that  the  requisite  integrals  must  exist,  and  that  the  singularities  in  the 
integrands  must  be  integrable,  it  is  instructive  to  demonstrate  this  afresh 
taking  into  account  the  explicit  knowledge  of  the  logarithmic  singularities. 
We  accordingly  digress  to  give  a  brief  proof  that  the  contribution  of  the 
singularity  of  K(m)  to  these  integrals  is  finite. 


When  m  is  very  close  to  unity,  nei,  the  elliptic  integral  E(m)  is  also 
nearly  unity  (and  therefore  finite),  but  the  elliptic  integral  of  the  first 
kind  has  the  behavior 


K(ra)  *  2  ,n(r^i)  (74) 

For  the  integrations  involved  in  the  expression  (56)  for  the  matrix  coef¬ 
ficient  ajn,  the  region  in  the  integration  plane  surrounding  the  singularity 
can  be  taken  to  be  the  strip  along  the  diagonal  described  by  r-A£s<r+Ar  where 
A  is  very  small.  In  this  region,  the  dominant  contributions  to  the  Cj(rls) 
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factors  in  the  integrand  are  identified  from  Eqs.  (64)  and  (69)  to  be 


C1(rls)  =  2rmK(m) ;  C2(rls)  =  4rK(m) 

If  one  assumes  that  the  derivatives  $'(r)  of  the  basis  functions  are  finite, 
then  the  corresponding  contributions  from  such  singular  terms  to  the  integrals 
in  Eq.  (56),  which  define  the  ajn,  can  be  generically  written 


{Contribution  to  a.  } 

J  ^ 


1  r+A 

A  A 

F(r,s)  r  K(m)  ds  dr 
0  r-A 


+  0(A) 


(75) 


where,  with  the  assumption  previously  stated,  F(r,s)  denotes  a  factor  that  is 
finite  in  the  domain  of  integration. 


Because  of  the  smallness  of  A,  the  mean  value  theorem  of  the  integral 
calculus  allows  us  to  neglect  fluctuations  in  F(r,s)  resulting  from  changing 
s,  so  we  replace  F(r,s)  by  F(r,r).  Then,  in  order  to  perform  the  integral  over 
the  domain  of  s,  we  change  the  variable  of  integration  to  £,  such  that 

2 

s=r(l+f),  l-m  «(§=£)  £2;  K(m)  ~  \  In  (64/^)  (76) 

Such  substitutions  transform  Eq.  (75)  to 
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1  A/r 


{Contribution  to  a.}  2  F(r,r)(2r2)  ln(8/f)  dfdr  +  0(A) 


(77a) 


0  0 


1 


F(r,r)  (2rA)ln(8r/A)  dr  +  0(A) 


(77b) 


0 


The  integrand  factor  (2rA)ln(8r/A)  is  bounded  in  absolute  value  for  all  com¬ 
binations  of  r  and  A,  given  that  both  are  between  0  and  1.  This  is  clear 
because  x[ln(l/x)]  is  less  than  e"^0.37.  Thus,  providing  F(r,r)  is  finite 
for  all  r  between  0  and  1,  the  integral  in  Eq.  (77b)  is  bounded,  and  it  goes 
to  zero  when  A  goes  to  zero. 

As  is  explained  further  below,  it  is  advantageous  to  use  basis  functions 
$(r)  that  go  to  zero  like  (l-r)l/2  when  r*l.  Such  would  require  $■'  to  have  an 
integrate  singularity  at  r=l.  However,  the  above  demonstration  of  bounded¬ 
ness  would  no  longer  be  valid  because  F(r,s)  would  be  singular  when  either  r 
or  s  are  unity.  Thus,  there  would  be  a  confluence  of  two  or  three  singular¬ 
ities  at  the  point  r=s=l.  To  check  whether  this  precludes  the  guaranteed 
existence  of  finite  matrix  coefficients,  it  is  sufficient  to  examine  the 
existence  of  the  integral 


1  1 


{Contribution  to  a.  }  2  M 

jnJ 


(78) 


1-6  1-6 


where  M  is  finite  and  e  is  small  compared  to  unity.  The  above  is  an  integral 
over  a  small  square  of  side  length  e  with  a  corner  at  the  point  r=s=l.  Alter¬ 
natively,  it  is  sufficient  to  consider  integration  over  the  quadrant  of  a 
circle  centered  at  r=s=l  and  having  radius  e.  With  the  latter  understanding, 
we  let  1-r  be  ae  cos^  and  1-s  be  ae  sin^,  such  that  ds  dr  is  replaced  by 
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e2a  da  and  a  ranges  from  0  to  1,  while  <f>  ranges  from  0  to  r/2.  The  quad¬ 
rant  integral  corresponding  to  the  examination  exercise  posed  by  Eq.  (78) 
above  then  becomes 


t/2  1 


{Contribution  to  a.  } 
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=  Me 


lnfevreVfcos^-stn^2!!  da  ^ 
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(79) 


Both  of  the  indicated  integrations  over  <f>  in  the  latter  expression  exist  and 
have  values  of  the  order  of  unity,  even  though  the  integrands  have  one  and  two 
singularities,  respectively.  Both  terms  are  finite  and,  moreover,  they  go  to 
zero  when  e  goes  to  zero. 


One  may  conclude  from  the  analysis  of  the  present  section  that  the  inte¬ 
grated  Green's  functions  may  be  evaluated  numerically  without  much  effort, 
even  though  the  original  integrals  defining  these  quantities  had  integrands 
with  singularities.  Both  of  these  integrated  Green's  functions  have  a 
logarithmic  singularity  at  r=s,  so  when  these  functions  appear  as  factors  in 
an  integrand,  as  in  Eq.  (56),  the  integration  scheme  must  not  explicitly  ask 
for  values  at  such  arguments.  However,  the  presence  of  these  singularities 
does  not  cause  any  of  the  anticipated  integrands  to  be  nonintegrable. 
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10.  EVALUATION  OF  THE  MATRIX  ELEMENTS 


Near  the  edge  of  the  disk  the  exact  solution  for  the  acoustic  pressure  on 
the  surface  must  behave  to  leading  order  in  the  distance  from  the  edge  in  the 
same  manner  as  would  a  solution  of  Laplace's  equation  near  a  knife  edge.  Thus 
if  a  is  radial  distance  from  the  edge  and  if  ^  is  angle  about  the  edge,  such 
that  f  is  0  on  the  front  side  and  2r  on  the  back  side,  then 

92p/9a2  +  a”*9p/9a  +  a~292p/9^2  ~  0  (80) 

for  sufficiently  small  a.  The  pressure  p  must  be  finite  near  the  edge,  and  it 
must  satisfy  the  rigid  wall  boundary  conditions 


dp/da  =0  at  y=0'and  at  ^=2t  (8i) 

This  "boundary  value  problem"  can  be  solved  by  the  method  of  separation  of 
variables,  with  the  result  that,  to  leading  order  in  a, 


P  -  A  +  B  o1/2cos(^/2)  (82) 

where  A  and  B  are  "constants".  For  the  particular  case  of  sound  generated  by 
a  disk  in  rigid  body  transverse  vibration,  symmetry  and  continuity  require 
that  p  be  zero  at  the  edge,  so  the  "constant"  A  is  identically  zero.  Hence, 
on  the  front  surface,  where  jf=0,  the  implication  of  the  above  equation  is 

p  ~  (distance  from  edge)1^2  (83) 
near  the  edge  of  the  disk. 

The  above  deduced  behavior  must  be  exhibited  by  the  exact  solution  for  p 
of  either  the  differential  integral  equation  or  of  the  variational  principle. 
It  is  not,  however,  a  requirement  that  must  be  imposed  at  the  outset  on  the 
trial  functions  or  the  basis  functions.  The  derived  approximate  solution,  if 
the  number  of  basis  functions  is  sufficiently  large,  can  be  expected,  in  the 
aggregate,  to  approximate  the  behavior  of  Eq.  (83). 
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Since  we  know  the  result  (83)  in  advance  of  a  detailed  numerical  solution, 
a  strong  argument  can  be  made  that  it  should  be  incorporated  into  the  trial 
functions  at  the  outset,  for  then  a  faster  convergence  toward  the  exact 
solution  might  be  achieved.  One  possibility  for  doing  this  is  to  include  a 
factor  (l-r2) 1/2  in  each  of  the  *n(r).  One  could  take,  for  example, 
$l(r)=(l-r2)l/2,  which  turns  out  to  be  the  exact  result  [Pierce,  1981]  for  the 
pressure  distribution  on  the  disk  in  the  limit  ka+O.  Doing  so,  however, 
introduces  another  singularity  into  one  of  the  integrands  of  the  integrals 
that  define  the  matrix  elements  ajn  in  Eq.  (56).  Clearly,  $i'(r)  is  infinite 
at  r=l.  The  singularity  is  integrable,  but  integration  over  singular  inte¬ 
grands  is  an  inherent  source  of  numerical  difficulties.  Such  difficulties  are 
often  surmountable  with  a  change  of  integration  variable  that  transforms  the 
original  integrand  into  one  in  which  the  singularity  does  not  appear. 

With  the  purpose  just  described  in  mind,  we  denote  the  desired  transfor¬ 
mation  as 


r  =  g(u);  s  =  g(w) 


Then 


(84) 


Wr>  dr  ‘  du  (85) 

If  g(u)  is  chosen  appropriately,  then  (d/du)$n(g(u))  will  be  finite  at  the 
value  of  u  where  g(u)=l,  even  if  %(l)  is  not.  It  is  convenient  to  have  the 
limits  of  the  domain  of  u  and  w  match  those  of  r  and  s,  so  we  introduce  the 
requirement  that 


9(0)  =  0;  g(l)  =  1  (86) 

The  choice  for  the  integration  variable  transformation  function  g  could  be 
different  for  different  basis  functions,  but  it  is  advantageous  for  us  to 
require  that  it  be  the  same  for  all  elements  within  the  same  set,  such  that  no 

symmetry  properties  are  lost;  this  is  assumed  to  have  been  done  in  the  discus¬ 
sion  that  follows. 
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In  the  case  of  the  basis  function  $i=(l-r2) 1/2,  we  shall  set 


g(u)  =  sin  [p]  (87) 

This  satisfies  Eqs*  (86),  and  the  derivative  in  Eq.  (85)  transforms  to 

^g(«)]  ■  U^1n2(f)]1/2 

3  '  f  s1"  (f)  (88) 

which  is  well-behaved  for  all  u. 

To  simplify  the  notation,  we  let  rn(u)  denote  the  n-th  basis  function  when 
it  is  regarded  as  a  function  of  the  alternate  integration  variable  u,  so  that 

rn(u)  =  ¥n[g(u)]  (89) 

Then  substitution  of  Eqs.  (84),  (85),  and  (89)  into  Eq.  (56)  yields 


1  1 

(ka)2  [ 

/  - 
0  0 


g' (u)g' (w)rj(u)rn(w)C1(g(u) lg(w))  dwdu 


1  1 

rj(a)r'(w)Ci(g(u)lg(w))  dwdu  (go) 

J  J 
0  0 

Because  the  source  variable  r  and  field  variable  s  undergo  the  same  trans¬ 
formation,  the  singularity  contributed  by  the  integrated  Green* s  functions 
continues  to  lie  on  the  diagonal,  which  is  now  the  line  w=u.  This  line  can  be 

taken  to  be  a  boundary  line  for  the  integration  domain,  by  using  the  simple 
identity 
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1  a 


1  1 
f  [ 

F(u,w)  dwdu 

J  4 
0  0 


[f(u,w)  +  F(w,u)jdw  du 


0  0 


(91) 


where  the  function  F(u,w)  is  arbitrary.  With  the  recognition  that  the  inte¬ 
grated  Green's  functions  are  unchanged  by  interchange  of  their  arguments,  the 
application  of  this  identity  to  Eq.  (90)  yields 


1  u 

[rj  (U)rn(w)+I,j  (w)rn(u)]g'  (u)g‘  (w)  C|(g(u)  lg(w))dw  du 

0  0 

1  u 

[rj (u)Fn (w)+rj (w)r^ (u) J  C2(g(u) lg(w))dw  du  (92) 

J  J 
0  0 

In  a  similar  manner,  the  "forcing  vector"  components  bn,  given  previously 

by  Eq.  (57),  when  expressed  in  terms  of  the  new  variable  of  integration  u, 
become 


1 

bn  =  2xika  f  g(u)g' (u)Tn(u)  du  (93) 

0 

The  integrations  required  to  evaluate  the  components  bn  may  be  effected  by 
direct  numerical  quadrature;  Simpson's  1/3  rule  should  be  adequate.  The 
interval  0<u<l  is  regarded  as  composed  of  J  intervals,  where  J  is  even;  then 

J 

bn  =  2xi ka>~  rn(uj.)g(uj.)g'  (if.)*.  (94) 

j=0 
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where  Uj  -  j/J  is  the  endpoint  of  the  j-th  integration  interval.  The  weighting 
coefficients  i/j  for  Simpson's  Hl/3"  rule  are  1/(3J),  4/(3J)r  2/(3J),  4/(3J), 

. .  4/ (3J) ,  and  1/(3J)  for  j  =0,  1,  2,  3, . .  J-l,  and  J.  The  J-th 

term  in  the  above  sum  is  identically  zero  because,  as  remarked  previously,  one 
must  have  rn(l)  =  0.  Recall  that  all  admissible  basis  functions  must  vanish 
at  the  edge  of  the  disk. 

The  double  integrations  required  for  the  evaluation  of  the  matrix  coef¬ 
ficients  ajn  in  Eq.  (90)  cannot  be  evaluated  by  sequential  application  of  a 
one-dimensional  integration  rule.  Such  an  approach,  which  would  discretize  u 
and  w  in  the  same  manner,  would  tacitly  assume  that  the  integrand  be  finite 
along  the  line  w=u  (where,  in  fact,  the  integrand  is  singular),  since  the 
integration  algorithm  would  use  values  of  the  integrand  evaluated  at  points 
along  this  line.  We  circumvent  this  difficulty  by  using  the  alternative 
expression,  Eq.  (92),  and  by  employing  an  area  integration  rule  based  on 
interior  points.  As  shown  in  Fig.  4,  the  triangular  integration  domain,  which 
is  bounded  by  the  three  lines  w=0,  u=l,  and  w=u,  is  segmented  into  K(K-l) 
slanted  squares,  and  2K  triangles,  based  on  K  divisions  along  the  lines  w=0 
and  u=l.  The  value  of  K  need  not  be  the  same  as  the  number  of  intervals  J 
that  are  used  in  the  evaluation  of  the  bn.  For  integration  over  each  of  the 
interior  squares,  one  uses  a  nine-point  integration  rule  [Davis  and  Polonsky, 
1972],  such  that  the  integral  over  a  square  is  approximated  by  the  area  of  the 
square  times  a  weighted  average  of  the  values  of  the  integrand  at  nine 
interior  points.  The  locations  of  these  points  and  the  weighting  factors  are 
selected  such  that  the  integration  will  be  exact  if  the  integrand  is  a  poly¬ 
nomial  of  the  fifth  degree.  The  relative  location  of  these  interior  points 
within  a  single  square  is  shown  in  Fig.  5. 

Let  2A  -  1/K  be  the  width  of  a  subdivision  interval  along  the  line  w=0 
between  u=0  and  u=l.  Each  of  the  interior  squares  will  then  have  length  2*/2a 
on  a  side,  diagonals  of  length  2A,  and  areas  2A2.  The  centerpoints  of  these 
interior  squares  are  at 
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w 


Figure  4.  Discretization  of  integration  domain  into  segments. 
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=  2kA  +  (m-l)A;  w, 


=  mA 


m  =  1 


I  •  •  •  f 


2K-2k 


ukm,c 


km,c 


k  =  1. 


•  •  •  I 


K-l 


(95) 


such  that  squares  distinguished  by  the  same  value  of  the  index  k  have  center- 
points  that  lie  along  a  common  line  parallel  to  the  diagonal  line  w=u;  squares 
with  the  same  k  are  numbered  by  the  index  m  in  the  order  in  which  they  are 
encountered  as  one  moves  up  this  line  toward  increasing  u  and  w.  The  coor¬ 
dinates  of  the  nine  integrand  sampling  points  within  a  given  square  and  the 
corresponding  weights  (which  sum  to  unity)  are 


■W")  -  Vc+  a(q)4: 

Wkm^  =  wkm,c+  b^A-  q  =  1.  2,  ...  9 

a(l)  =  0 

b(l)  =  0 

a(2)  =  -  (3/20) 1/2 

b(2)  =  +  (3/20) 1/2 

a(3)  =  +  (3/20) 1/2 

b(3)  =  +  (3/20) 1/2  . 

a(4)  =  +  (3/20) 1/2 

b(4)  =  -  (3/20) 1/2 

a(5)  =  -  (3/20) 1/2 

b(5)  =  -  (3/20) 1/2 

a(6)  =  -  (3/5) 1/2 

b(6)  =  0 

O 

II 

b(7)  =  +  (3/5) 1/2 

a (8)  =  +  (3/5) 1/2 

b  (8)  =  0 

O 

II 

o> 

b(9)  =  -  (3/5) 1/2 

"sqd)  *  16/81:  "sq(2) 

"sqd)  -  "sqd)  ‘  "sq® 

=  I/sq(3)  =  I/sq(4)  =  ^sq^5  =  10/81 
"  "sqW  =  25/324  (96) 

The  integrations  over  the  K  triangles  that  fill  the  gaps  left  by  the 
squares  along  the  line  w=0  (lower  edge,  abbreviated  le)  can  be  similarly  per¬ 
formed  using  a  seven-point  integration  rule.  The  lower  edges  of  these 
triangles  have  their  centers  at  u  =  A,  3A,  5A,  (2K-1)A;  each  has  height  A 
and  area  A^.  The  coordinates  of  the  integrand  sampling  points  within  the  k-th 

triangle  and  the  corresponding  weights  are  derived  by  mapping  the  45°-45°-90° 
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triangle  to  an  equilateral  triangle  and  then  using  a  standard  integration 
scheme  for  the  latter.  [One  of  the  expressions  on  page  483  in  the  Dover 
edition  of  Abramowitz  and  Stegun's  Handbook  of  Mathematical  Functions  has  an 
error.  What  appears  there  as  -(15)1/2+1  should  be  -[(15)1/2+1] .]jhe  result  of 
such  a  derivation  yields 

ule,k^  =  (2k_1)A  +  a^A;  wle,k(q)  =  £(q)A'  Q  3  1»  2,  ..,  7 


«(1) 

=  0 

P(  1) 

= 

1/3 

a(  2) 

=  0 

P(2) 

= 

[9  + 

2(15) 1/2]/21 

a(3) 

=  - 

[1  + 

(15) 

1/2]/7 

P(  3) 

= 

[6  - 

(15)1/2]/21 

a(4) 

=  + 

[1  + 

(15) 

1/2]/7 

P(  4) 

= 

[6  - 

(15) 1/2]/21 

a(5) 

=  0 

P(  5) 

S 

[9  - 

2(15) 1/2]/21 

a(  6) 

=  - 

[(15) 

1/2. 

■13/7 

P(  6) 

= 

[6  + 

(15) 1/2]/21 

a(7) 

=  + 

[(15) 

1/2. 

13/7 

;(7) 

= 

[6  + 

(15) 1/2]/21 

ytri 

(1)  = 

270 

1200 

I 

II 

c\T 

u 

4-> 

^tri 

(3) 

=  ~  . (4)  =  155.-  (15) 

tri v  '  1200 

^tri 

(5)  = 

vtr1 

(6) 

*  "trl <7> 

.  155  +  (15) 
1200 

1/2 

Similarly,  for  the  K  triangles  with  hypotenuses  along  the  line  u=l  (right 
edge,  abbreviated  re),  one  has 

ure,k^  =  1  ~  ^(q)A;  wre,k(q)  =  (2k“l)A  +  ®(q)A;  q  =  1,  2,  ..,  7  (98) 

with  the  a(q),  p(q) ,  and  i/-tri  (q)  the  same  as  in  Eq.  (97). 

The  integration  scheme  just  described  gives  the  following  approximate 
representation  for  the  double  integral  of  any  function  F(u,w)  over  the  unit 
square: 
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1  1 
A  A 

F(u,w)  dwdu 

0  0 


1  u 

A  A 

[f(u,w)  +  F(w,u)Jdw  du 

J  4 
0  0 


K-l  2K-2k  9 

■  Z  Z  X2A2s(ukn,W-  «la,(D)  "sqW 

k=l  m=l  q=l 
K  7 

+  Z 

k=l  q=l 
K  7 

+  Z  (") 

k=l  q=l 

where 

A  =  1/(2K);  S(utw)  =  F(u,w)  +  F(w,u)  (100) 

For  the  computation  of  a  particular  ajn,  the  function  F(u,w)  is  identified  as 


Integrand  F(u,w)  for  ajn=  (ka)2^.  (u)rn(w)g'  (u)g1  (w)  C1(g(u)  I g (w) ) 

"  rj Cu)rn (w)  c2(g(u) ig(w))  (ioi) 


Many  numerical  operations  are  required  to  evaluate  the  matrix  elements  ajn 
using  Eqs.  (99-101)  for  all  combinations  of  j  and  n.  One  saving  comes  from 
the  recognition  that  the  sampled  values  of  the  integrated  Green's  functions 
appearing  there  are  independent  of  the  (j,n)  pair  under  consideration.  Hence, 
the  values  of  g1 (u)g' (w)Ci(g(u) lg(w))  and  C2 (g (u) lg(w))  for  all  points  (u,w) 
used  in  the  sums  in  Eq.  (99)  may  be  stored  in  arrays,  and  then  recalled  when¬ 
ever  necessary  to  calculate  the  integrand  (101)  at  the  corresponding  points. 
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11.  FINITE  ELEMENT  BASIS  FUNCTIONS 


The  implication  thus  far  has  been  that  the  set  of  basis  functions  consists 
of  a  few  analytical  functions  satisfying  the  boundary  condition  p(l)=0. 
However,  the  generic  trial  functions  will  also  be  admissible  if  the  basis 
functions  are  not  smooth  curves,  but  merely  continuous,  with  a  finite  number 
of  derivative  discontinuities. 

Suppose  one  seeks  to  represent  the  surface  pressure  distribution  p(r)  by  N 
discrete  values  at  evenly  spaced  radial  distances  rn=(n-l)/N.  As  shown  in 
Fig.  6,  a  linear  interpolation  between  such  discrete  values  may  be  represented 
as  a  superposition  of  piecewise  linear  functions.  Let  $n(r)  denote  a  linear 
finite  element  basis  function  centered  at  rn.  We  define  this  basis  function 
to  be  unity  at  r=rj,  and  to  vanish  beyond  the  adjacent  points  rn_i  and  rn+i. 
A  simple  algorithm  for  evaluating  these  basis  functions  defines  the  local 
di stance 


d=r-rn=r-  (n-l)A;  A  =  1/N 


and  then  sets 


(102) 


fl-ldl/A 

lo 


if  I d I <A 
if  1  d 1 >A 


»„(r) 


’-1/A 

1/A 

.0 


if  0Sd<A 
if  -A<d<0 
if  I d I >A 


(103) 


When  we  use  these  finite  element  basis  functions  to  form  a  generic  trial  func¬ 
tion  as  in  Eq.  (54),  it  is  immediately  apparent  that  the  coeficients  of  the 
basis  functions  are  identical  to  the  discrete  pressure  values,  that  is, 
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Figure  6a.  Linear  interpolation  of  a  discretized  pressure. 


Figure  6b.  Finite  element  basis  functions. 
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(104) 


A 

Pn  =  P^rn^  where  rn  =  (n_1)A 

In  principle,  we  could  employ  the  finite  element  basis  functions  directly 
using  the  general  formulation  and  numerical  scheme  developed  in  the  preceding 
section.  However,  doing  so  would  lose  potential  computational  savings  result¬ 
ing  from  the  fact  that  the  finite  element  basis  functions  vanish  over  a  large 
portion  of  the  range  0£r£l.  Such  savings  may  be  very  significant,  because 
many  basis  functions  may  be  required  for  an  accurate  description  of  the  pres¬ 
sure  distribution. 

For  such  a  purpose,  it  is  first  convenient  to  split  the  overall  sum  of  the 
type  in  Eq.  (99)  for  the  evaluation  of  a  given  ajn  using  Eqs.  (100)  and  (101) 
into  two  partial  sums  Sjn  and  Sn j ,  such  that 


where 


ajn  =  anj  Sjn  +  Snj 


Sjn  ‘ 


K-l  2K-2k  9 

k=l  m=l  q=l 
K  7 

+  Z  Z  *2  FJn(ule,k(q)-  wle,k(q))  ytH  W 

k=l  q=l 
K  7 

+  Z  ZA2Fon<V,kW. 

k=l  q=l 


Fjn(u,w)  =  (ka^iyu^Mg-Mg-tw)  C1(g(u)  lg(w)) 

-  rj (u)r^(w)  C2(g(u)lg(w)) 


(105) 


(106) 


(107) 
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Since  the  sum  in  Eq.  (106)  uses  only  points  for  which  the  w  coordinate  is 
less  than  the  u  coordinate,  the  sum  defining  the  Sjn  will  be  identically  zero 
if  all  of  the  u's  for  which  Tj(u)  is  nonzero  are  less  than  all  of  the  u's  for 
which  rn(u)  is  nonzero. 


By  definition,  Tj(u)  =  ¥j(g(u)).  For  the  finite  element  basis  functions 
defined  by  Eqs.  (103),  the  products  $j(r)$n(s)  and  $j(r)$^,(s)  are  nonzero  only 
in  the  domain  Rjn  defined  by 


0  S 


111 

N 


n  <r  n-2  .  e  <  n 
0  ^  S  s  S  N 


(108) 


This  domain  is  depicted  in  Fig.  7a.  Because  the  boundaries  of  Rjp  correspond 
to  constant  values  of  r  or  s,  the  transformation  of  r=g(u)  and  s=g(w)  maps  rjn 

into  the  region  Ujn  in  the  u-w  plane,  as  depicted  in  Fig.  7b.  This  region  is 
defi ned  by 


0  S 


i  u  S 


where  g~l  denotes  the  inverse  transformation. 


(109) 


Figure  7b  also  depicts  the  region  covered  by  several  integration  area- 
elements  (i.e.,  those  squares  and  triangles  used  in  the  breaking  up  of  the 
overall  integration  domain  in  the  u-w  plane).  Any  area-element  that  does  not 
overlap  the  finite  element  domain  Ujp  may  be  skipped  in  the  formation  of  the 
partial  sum  Sjn.  Furthermore,  when  an  area  element  only  partially  overlaps 
Ujn»  The  contributions  from  integrands  evaluated  at  points  (u,w)  that  do  not 
satisfy  Eq.  (109)  may  also  be  skipped.  This  observation  reduces  the  number  of 
computations  because  it  avoids  operations  on  terms  that  would  eventually  be 
found  to  vanish. 
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s 


Figure  7a.  Domain  of  a  finite  element  in  physical  variables. 


Figure  7b.  Domain  of  a  finite  element  in  integration  variables. 
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Because  Eqs.  (105),  (106),  and  (107)  are  equivalent  to  Eqs.  (99),  (100), 
and  (101),  we  may  also  employ  the  partial  sums  to  evaluate  the  coefficients 
ajn  corresponding  to  distributed  basis  functions.  Hence,  the  only  difference 
in  the  algorithms  for  distributed  and  finite  element  basis  functions  is  that 
the  latter  should  be  checked  for  overlap  according  to  Eq.  (109).  We  also 
implemented  an  analogous  overlap  check  in  the  computation  of  the  summation  for 
bn,  Eq.  (94),  in  the  case  of  finite  element  basis  elements. 

12.  NUMERICAL  RESULTS 

An  analytical  solution  for  the  pressure  distribution  along  the  face  of  a 
thin  disk  is  available  for  the  limiting  case  of  very  low  frequencies,  ka+0. 
Thus,  we  chose  this  case  for  the  first  validation  of  the  variational  principle 
and  its  implementation.  When  the  pressure  on  the  positive  z  face  is  v0 
exp(-iwt) ,  with  ka  «  1,  then  the  pressure  is 

p  =  "  r  1>cv0ka(1  •  r2)1/2  (110) 

In  addition  to  providing  a  simple  expression  for  comparison  with  numerical 
results,  this  form  provides  an  indirect  method  of  verification.  If  this  mode 
is  used  as  one  of  several  modes,  the  amplitudes  of  the  others  should  be  very 
small  when  ka  «  1.  We  tested  the  accuracy  of  this  hypothesis  by  redefining 
the  modal  amplitudes  to  be  Pj /ka  and  selecting  the  modes  to  be 

=  r 2^J  "  ^(1  -  r2)  1/2  (HI) 

Table  1  gives  the  modal  amplitudes  as  a  function  of  the  number  of  modes  N  in 
the  series  expansion.  These  results  were  developed  by  applying  the  sine 
transformation  for  the  radial  distance,  Eq.  (87),  with  the  number  of  integra¬ 
tion  segments  set  at  K  =  10.  The  accuracy  of  the  variational  principle  is 
remarkable.  Note  that  the  first  mode  is  very  close  to  the  value  in  Eq.  (110), 
while  the  other  modes  are  three  orders  of  magnitude  smaller. 
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Table  1 


Approximate 

solutions  at 

ka  ->  0 

Coefficients  corresponding  to  the  shape  functions 

2(n-l)  i - 

r  ’  1-r2  ,  n  =  1,  2,  ...  N 

C1 

C2 

C3 

N  =  1 

-0.63875 

N  =  2 

-0.63875 

-6 

0.27369  X  10 

N  =  3 

-0.63875 

-0.21376  X  10"5 

0.27127  X  10-5 

Analytical 
solution  at 
ka  ->  0 

-0.63662 

Rayleigh-Ritz  approximations  of  the  acoustic  pressure  on  the  surface 
of  a  transversely  vibrating  rigid  circular  disk.  The  number  of  inte¬ 
gration  intervals  is  set  at  20  in  each  case. 
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We  also  checked  the  finite  element  results  for  ka+0  against  Eq.  (110). 
For  this,  we  must  consider  the  general  effect  of  inadequate  choices  for  the 
number  of  elements  N  and  the  number  of  integration  segments  K.  The  value  of  N 
must  be  sufficiently  large  to  represent  the  pressure  as  a  sequence  of  discrete 
values.  We  estimated  that  N  £  5  was  appropriate  for  the  pressure  given  by 
Eq.  (110).  When  the  number  of  elements  has  been  set,  the  value  of  K  must  be 
sufficient  to  place  several  integration  points  within  the  domain  covered  by  a 
finite  element.  We  found  that  sensible  results  are  obtained  if  K  >  N. 

The  transformation  in  Eq.  (87)  was  used  in  the  evaluation  of  the  finite 
elements  model,  although  it  is  not  necessary  to  do  so.  We  will  discuss  the 
reasons  for  this  choice  later.  Comparison  of  finite  element  predictions  for 
several  values  of  K  when  N  =  5  and  10  are  shown  in  Figure  8  and  9,  respec¬ 
tively.  The  agreement  between  the  finite  element  formulation  and  the 
analytical  solution  is  evident.  Clearly,  ten  elements  provide  better  resolu¬ 
tion  of  the  pressure  distribution,  but  the  results  exhibit  substantial 
numerical  noise  with  increasing  K. 

Such  noise  could  have  two  causes.  The  computations  were  performed  on  a 
VAX  11/750,  which  due  to  hardware  limitations  in  its  current  configuration 
slows  drastically  when  performing  double  precision  complex  arithmetic  in 
programs  that  access  virtual  memory.  However,  we  checked  the  results  with 
comparable  calculations  on  a  CDC  CYBER  785,  which  doubles  the  arithmetic 
precision  because  the  word  size  on  the  CYBER  is  twice  as  large  as  that  on  the 
VAX.  Nevertheless,  we  found  that  the  results  changed  very  little.  We  have 
concluded  that  the  primary  source  of  error  generation  for  finite  elements  is 
the  usage  of  the  piecewise  linear  modes  in  Eqs.  (103).  As  the  interval 
covered  by  each  element  decreases,  continuous  analytical  derivatives  are 
modeled  in  those  modes  as  large  values  that  change  sign  at  the  center  point. 
When  ka  =  0,  the  integrand  forming  the  coefficient  array  [A]  depends  solely  on 
these  derivatives,  so  it  is  not  surprising  to  see  the  predictions  display 
numerical  noise.  We  believe  that  these  problems  would  be  ameliorated  through 
the  application  of  smoother  finite  elements,  such  as  polynomials. 
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o 


Figure  8.  Pressure  distribution  when  ka=0  using  five  finite  elements 
(N=5) . - :  K=5 ;  . :  K=10; - :  K=20;  -  :  Eq.  (110). 
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o 


Figure  9.  Pressure  distribution  when  ka=0  using  ten  finite  elements  (N=10). 
- :  K=10; - :  K=20;  - :  Eq.  (110). 
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The  modes  in  Eq.  (Ill)  are  such  that  one  matches  the  analytical  solution. 
A  fairer  test  is  the  ability  of  the  variational  principle  to  predict  the 
correct  pressure  distribution  when  the  modes  that  are  selected  merely  satisfy 
the  boundary  condition  that  p  =  0  at  the  edge.  Our  choice  for  this  evaluation 
was  a  half-range  cosine  series, 


*j a  cos  [((2J-i)r]  (112) 

Figure  10  shows  that  an  expansion  using  five  of  these  modes  compares  quite 

favorably  with  Eq.  (110).  Indeed,  the  prediction  is  virtually  identical  to  a 

Fourier  series  expansion  of  Eq.  (110)  using  Eq.  (112). 

The  results  in  Figure  10  were  computed  using  the  transformation  in 

Eq.  (87).  Such  a  transformation  is  not  required  for  these  functions,  nor  for 
the  finite  element  modes.  Both  Eq.  (112)  and  Eq.  (103)  give  finite  values  for 
at  r=l.  Hence  it  would  be  possible  in  either  case  to  perform  the  numerical 
integrations  leading  to  [A]  with  an  identity  transformation  of  g(u)=l. 

When  we  employed  the  identity  transformation,  we  found  that  favorable 
comparison  with  the  analytical  solution  required  more  modes  or  finite  ele¬ 
ments.  This  observation  has  a  ready  explanation.  The  integration  scheme 
described  by  Figures  4  and  5  uses  a  uniform  mesh  in  the  u-w  plane.  The  trans¬ 
formation  in  Eq.  (87)  maps  these  points  into  a  nonuniform  mesh  in  the  r-s 

plane,  such  that  the  density  of  integration  points  increases  monotonically  as 
r+1  and  s+1.  Thus,  the  transformation  gives  a  better  description  of,  and 

greater  emphasis  on,  the  behavior  near  the  edge,  where  diffraction  effects  are 

most  significant. 

Based  on  these  observations,  we  decided  to  perform  all  evaluations  of  the 
pressure  distribution  for  non-zero  ka  using  Eq.  (87).  Also,  the  nature  of  the 
slope  singularity  at  r=l  is  not  altered  when  ka^O.  For  this  reason,  further 
evaluations  using  analytical  modes  shall  only  be  based  on  Eq.  (111). 
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Figure  10.  Pressure  distribution  when  ka=0  using  sinusoidal  nodes: 
- :  N=5  &  K=10;  - :  Eq.  (110). 
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Data  with  which  the  results  of  the  variational  principle  may  be  compared 
is  somewhat  sparse.  Fortunately,  the  problem  of  diffraction  of  sound  by  a 
circular  disk  Is  analogous  to  the  present  radiation  problem.  Consider  a  plane 
wave  at  normal  incidence  to  a  thin  circular  disk  that  is  fixed  in  space.  Let 
ps  denote  the  scattered  acoustic  wave,  and  let  the  z  axis  coincide  with  the 
axis  of  symmetry.  For  a  wave  incident  on  the  face  z=0+,  corresponding  to 
propagation  in  the  direction  of  decreasing  z,  the  total  acoustic  pressure  is 

P  =  [Ps(r,z)  +  p0exp(-ikaz)]exp  (-iwt)  (113) 

where  r  and  z  are  nondimensional  cylindrical  coordinates  and  p0  is  the  ampli¬ 
tude  of  the  incident  wave. 


Since  the  disk  is  stationary,  the  acceleration  of  a  fluid  particle  must 
vanish  at  z=o,  which  leads  to 


ika  n 


(114) 


For  comparison,  we  find  in  Eq.  (4)  that  the  boundary  condition  for  acoustic 
radiation  from  a  disk  is 


ai  I  z=0  =  1w'vo  (115) 

when  the  disk  velocity  is  v0  exp  (-iwt).  Comparison  of  Eqs.  (114)  and  (115) 
shows  that  results  for  scattering  may  be  applied  to  the  radiation  problem  by 
letting  p0  =  pcv0. 


The  diffraction  problem  was  solved  by  Leitner  [1949]  using  oblate 
spheroidal  wave  functions.  Results  from  that  analysis  were  obtained  for 
ka=l,2,3,4,  and  5,  so  those  are  the  cases  we  shall  consider  here.  When  ka+0, 
the  pressure  amplitude  is  a  negative  imaginary  quantity,  corresponding  to 
pressure  that  is  in  phase  with  the  acceleration,  in  other  words,  a  virtual 
mass  impedance.  When  ka  y  0,  the  pressure  has  both  real  and  imaginary  parts. 
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Figures  11  -  15  compare  the  results  for  the  analytical  modes  in  Eq.  (Ill) 
with  those  obtained  by  Leitner.  The  analytical  mode  results  were  obtained  for 
series  truncation  at  N=5,  while  10  finite  elements  were  used. 

The  analytical  modes  give  excellent  predictions.  Indeed,  the  slope 
discontinuity  at  the  edge  seems  to  be  modeled  better  by  the  variational 
principle.  Another  highlight  of  Figures  13  and  14  is  that  they  show  that  the 
earlier  results  for  ka  =  1  and  2  were  interchanged. 

The  finite  element  model  seems  to  work  reasonably  well.  However  the 
accuracy  is  not  as  good  as  that  obtained  from  the  analytical  modes.  One  might 
think  that  the  discrepancies  are  due  to  the  relatively  small  number  of 
elements  and  integration  points.  Figure  16  shows  the  pressure  distribution 
for  ka=5  using  more  elements  and  a  finer  integration  mesh.  Little  improvement 
is  obtained  by  increasing  the  number  of  integration  points,  whereas  increasing 
the  number  of  elements  increases  the  numerical  noise  associated  with  the 
piecewise  linear  nature  of  the  elements. 

The  numerous  computer  runs  we  made  led  us  to  inquire  about  the  relation¬ 
ship  between  the  execution  time  and  the  values  of  K  and  N.  We  derived  an 
estimate  from  the  recognition  that  the  bulk  of  the  computations  are  associated 
with  three  distinct  operations: 

(1)  evaluation  of  the  integrated  Green's  functions  Cj (gCu^i ) I g (w^t ) ) , 

(2)  formation  of  the  coefficients  ajn, 

(3)  solution  of  the  simultaneous  equations  to  find  the  modal  amplitudes. 

It  is  clear  from  Figure  4  that  the  number  of  integration  points  is  propor¬ 
tional  to  K2,  so  the  time  required  for  phase  1  may  be  expected  to  have  a 
comparable  dependency.  Note  that  the  number  of  modes  is  irrelevant  for  this 
task.  For  the  second  phase,  Eqs.  (106)  and  (107)  indicate  that  each  coef¬ 
ficient  ajn  is  obtained  from  a  summary  over  K2  terms.  Since  the  number  of 
these  coefficients  is  proportioned  to  N2,  we  expect  that  the  time  in  the 
second  phase  will  be  proportional  to  N2K2.  Finally,  the  time  required  to 
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Figure  11.  Pressure  distribution  when  ka=l. 

- :  analytical  modes  (N=5,  K=10); 

- :  finite  element  modes  (N=10,  K=20); 

. :  Leitner  [1949], 
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Figure  12.  Pressure  distribution  when  ka=2. 

- :  analytical  modes  (N=5,  K=10); 

- :  finite  element  modes  (N=10,  K=20); 

. :  Leitner  [1949], 
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Figure  13.  Pressure  distribution  when  ka=3. 

- :  analytical  modes  (N=5,  K=10); 

- :  finite  element  modes  (N=10,  K=20; 

. :  Leitner  [1949]. 
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Figure  14.  Pressure  distribution  when  ka=9. 

- :  analytical  modes  (N=5,  K=10); 

- :  finite  element  modes  (N=10,  K=20; 

. :  Leitner  [1949]. 
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Figure  15.  Pressure  distribution  when  ka=5. 

- :  analytical  modes  N=5,  K=10); 

- :  finite  element  modes  (N=10,  K=20; 

. :  Leitner  [1949]. 
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PRESSURE  AMPLITUDE 


Figure  16.  Influence  of  number  of  modes  and 
integration  points  using  finite  elements, 

ka=5. - :  N=10  &  K=20) ; - ;  N=10  &  K=40; 

. :  N=20  &  K=40. 
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solve  the  system  equations  is  independent  of  K.  The  complex  equation  solver 
we  employed  (LEQ2C  in  the  IMSL  library)  seems  to  require  N2  operation. 

These  heuristic  evaluations  led  us  to  expect  that  the  computation  time 
would  fit  T=CiK2+C2N2K2+C3N2.  Reasonably  close  fits  for  the  runs  we  made  on  a 
VAX  17/750  (with  a  floating  point  processor  but  no  array  processor)  are 

analytical  modes: 

T(sec)  st  0.75K2  +  0.01K2N2  +  0.14N2 
finite  elements 

T(sec)  a  0.66  K2  +  0.14  N2  (93) 

The  differences  between  these  estimates  are  attributable  to  the  reduction  in 
the  number  of  operations  introduced  by  the  overlap  check.  For  example,  the 
time  in  phase  (2)  is  not  dependent  on  K2N2  because  the  area  covered  by  a 
finite  element  pair  in  the  (r,s)  integration  domain  is  inversely  proportional 
to  N2. 
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8.  CONCLUSIONS 


We  have  derived  a  variational  principle  for  the  pressure  on  the  exterior 
surface  of  an  arbitrary  body  in  harmonic  motion.  The  principle  which  was 
obtained  from  the  Kirchhoff-Helmholtz  integral  theorem,  relates  the  surface 
pressure  directly  to  the  normal  component  of  velocity.  It  avoids  the  need  to 
solve  the  wave  equation  for  the  exterior  domain  subject  to  compatibility 
conditions  at  the  surface.  Current  approximate  methods  for  decoupling  the 
surface  pressure  from  the  exterior  domain,  such  as  the  doubly  asymptotic 
approximation,  rely  on  assumptions  regarding  the  relationship  between  pressure 
and  particle  velocity.  In  contrast,  the  variational  principle  relies  on  no 
assumptions  beyond  those  associated  with  linearized  acoustic  theory. 

We  specialized  the  variational  principle  to  the  case  of  a  thin  circular 
disk  vibrating  transversely.  Comparable  techniques  may  be  employed  to 
describe  other  thin  bodies,  such  as  curved  panels  and  bars.  We  developed  two 
techniques  for  applying  the  variational  principle  to  the  evaluation  of  the 
surface  pressure.  Both  involved  expansion  of  the  pressure  in  a  modal  series, 
with  the  difference  being  whether  the  modes  are  analytical  functions  covering 
the  entire  surface,  or  finite  elements  covering  only  a  segment  of  the  surface. 

The  results  obtained  by  using  analytical  modes  were  found  to  be  remarkably 
accurate.  Trial  functions  that  match  the  singularities  encountered  in  dif¬ 
fraction  of  an  edge  give  the  best  results,  but  reasonable  convergence  was 
obtained  for  more  general  functions.  The  finite  element  formulation  predicted 
the  pressure  less  accurately.  We  concluded  from  computations  for  a  variety  of 
parameters  that  piecewise  linear  finite  elements  generate  significant  numeri¬ 
cal  noise  when  the  number  of  elements  is  increased.  We  expect  that  the  usage 
of  polynomial  elements  would  significantly  improve  the  results. 

We  employed  both  formulations  to  obtain  the  pressure  for  nondimensional 
radii  in  the  range  kaS5.  We  could  have  presented  results  for  ka>5,  but  have 
not  done  so  because  of  a  lack  of  data  to  use  for  a  comparison.  The  trends 
suggested  by  the  present  results  are  that  higher  ka  would  require  more 
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analytical  modes  or  elements  in  order  to  accommodate  the  more  rapid  spatial 
variation  of  pressure.  Correspondingly,  the  mesh  of  integration  points  would 
have  to  be  finer. 

The  variational  principle  may  be  extended  to  other  geometries.  One  of  our 
current  efforts  is  devoted  to  a  finite  length  cylinder  in  axial  vibration. 
Systems  lacking  axi symmetry,  such  as  a  transversely  vibrating  cylinder,  may  be 
treated  by  these  techniques,  except  that  more  modes  and  more  numerical  inte¬ 
gration  points  would  be  required  to  treat  the  higher  dimensionality  of  such 
systems. 

Another  of  our  current  efforts  is  devoted  to  coupling  the  variational 
principle  to  a  modal  description  of  structural  vibration.  The  result  will  be 
a  highly  accurate,  yet  computationally  simple,  description  of  fluid-structure 
interaction.  We  have  also  shown  that  the  variational  principle  may  be  readily 
modified  to  describe  the  scattering  of  acoustic  waves. 
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